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1  Summary 

The  objective  of  this  project  was  to  develop  robust  numerical  methods  for  solv¬ 
ing  mathematical  models  of  nonlinear  phenomena  such  as  nonlinear  conservation 
laws,  surface/image/data  reconstruction  problems,  advection-dominated  flows,  multi¬ 
phase  flows,  and  free-boundary  problems,  where  shocks,  fronts,  and  contact  dis¬ 
continuities  are  driving  features  and  pose  significant  difficulties  for  traditional  nu¬ 
merical  methods. 

The  main  thrust  of  this  research  program  was  to  explain  some  intriguing  nu¬ 
merical  observations  reported  by  Lavery,  Jiang,  and  Guermond  [4] 1  that  seemed 
to  indicate  that  for  some  classes  of  PDE’s  equipped  with  non-smooth  coefficients 
and/or  non-smooth  right-hand  sides  it  pays  off  to  approximate  the  solution  directly 
in  L1.  Contrary  to  standard  stabilized  L2-based  techniques,  L1  -based  methods  did 
not  seem  to  require  additional  ad  hoc  tunable  coefficients  or  limiting  procedures. 

We  have  finally  elucidated  some  of  the  above  issues  and  we  can  report  that  the 
major  achievements  of  our  2.5  year  long  L1  -program  spanning  from  June  2009  to 
November  2011  are  the  following: 

(1)  L1  -minimization  for  PDEs:  We  have  shown  that  it  does  pay  off  to  work  in 
L 1  for  steady  equations  Hamilton- Jacobi  equations.  We  have  proved  in  par¬ 
ticular  that  L1  -based  finite  element  approximations  converge  to  the  physically 
relevant  solution,  i.e.,  the  so-called  viscosity  solutions.  We  have  developed  al¬ 
gorithms  for  solving  L1 -based  discrete  minimization  problems  and  we  have 

'Only  the  papers  written  by  the  Pis  under  the  umbrella  of  the  present  grant  are  cited  and  listed 
at  the  end  of  the  report  in  the  Reference  section. 
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proved  that  our  algorithms  converge  and  their  computational  complexity  is 
optimal  in  one  space  dimension. 

(2)  L1  -reconstruction  of  data:  We  have  developed  algorithms  for  L1  -based  data 
fitting  and  surface  reconstruction.  We  have  developed  a  localized  (easy  to 
parallelize)  version  of  the  L1-based  minimization  methods. 

(3)  Entropy  viscosity  for  nonlinear  conservation  laws:  As  a  by-product  of  the 
research  program  on  the  L1 -approximation  of  Hamilton  Jacobi  equations  we 
have  discovered  that  time-dependent  nonlinear  conservation  equations  can  be 
stabilized  by  using  the  so-called  entropy  viscosity. 

We  believe  that  the  entropy  viscosity  idea  is  an  important  conceptual  break¬ 
through.  The  main  stabilization  mechanism  in  this  method  is  a  nonlinear  dissipa¬ 
tion  proportional  to  the  local  size  of  an  entropy  production.  For  this  reason  we 
call  the  method  entropy  viscosity.  We  propose  to  use  the  local  residual  of  an  en¬ 
tropy  equation  to  construct  the  artificial  viscosity.  One  immediate  consequence  of 
this  choice  is  that  the  viscosity  is  proportional  to  the  entropy  production,  which  is 
known  to  be  large  in  shocks  and  to  be  zero  in  contact  discontinuities.  As  a  result, 
this  strategy  automatically  distinguishes  shocks  and  contact  discontinuities,  i.e., 
no  detection  of  contacts  and  no  artificial  compression  is  needed.  This  strategy  is 
inspired  from  the  observations  that  we  made  when  solving  Hamilton  Jacobi  equa¬ 
tions  using  the  L1  -minimization  technique:  approximate  L'-minimizers  solve  the 
PDE  in  the  smooth  regions  but  they  do  not  in  the  non-smooth  regions;  the  approx¬ 
imation  mechanism  is  dominated  completely  by  the  entropy  production  in  the 
non-smooth  regions.  We  have  implemented  this  idea  and  we  have  obtained  very 
encouraging  numerical  results.  In  particular,  we  have  observed  excellent  resolu¬ 
tion  of  contact  waves,  and  the  additional  cost  of  the  L1 -minimization  is  avoided. 


2  Hamilton-Jacobi  equations 

2.1  One  space  dimension 

We  have  introduced  an  L1  -based  minimization  technique  to  solve  steady  Hamilton- 
Jacobi  equations  in  [1 1],  We  have  proved  in  this  paper  that  the  method  converges 
to  the  unique  viscosity  solution.  Adding  an  entropy  penalty  to  the  L 1  functional 
is  critical  to  achieve  convergence  to  the  correct  solution. 

We  have  also  proposed  an  algorithm  of  optimal  complexity  to  solve  one¬ 
dimensional  steady  Hamilton-Jacobi  equations  in  [13]  and  [6].  We  have  proved 
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convergence  in  the  case  of  one-dimensional  stationary  Hamilton- Jacobi  equations 
with  convex  Hamiltonian.  Actually  we  have  proved  in  one  space  dimension  that 
our  algorithm  is  quite  similar  to  the  well-known  Fast  marching/Fast  sweeping 
methods  of  Sethian  and  Osher  which  means  that  the  Fast  marching/Fast  sweeping 
methods  construct  almost  minimizers  in  L1.  This  property  was  totally  unknown. 
We  conjecture  that  this  analogy  still  holds  in  higher  space  dimensions.  If  true  our 
methodology  would  then  prove  convergence  of  the  Fast  marching/Fast  sweeping 
algorithms,  since  we  have  already  proved  that  sequences  of  almost  L1  -minimizers 
converge  to  the  viscosity  solution  of  the  equation. 

We  show  in  Figure  1  the  solution  of 

j^\u\x)\  —  |  cos(27ra;)|  =  0,  u(0)  =  0,  u{l)  =  0,  (1) 

using  the  algorithm  developed  in  [13]  and  [6]. 


Figure  1:  Solution  to  (1)  using  100  degrees  of  freedom. 

2.2  Higher-space  dimensions 

We  have  revisited  the  one-dimensional  theory  developed  in  [11]  and  extended  it 
to  higher-space  dimensions  in  [12],  We  introduce  an  L1 -functional  augmented 
with  an  entropy  functional  that  is  dimension-dependent.  We  have  proved  in  [12] 
that  the  approximate  solution  obtained  by  minimizing  this  combined  functional 
converges  to  the  unique  viscosity  solution. 

The  method  has  been  tested  numerically.  We  show  in  Figure  2  two  compu¬ 
tations  performed  on  the  eikonal  equation  using  two  types  of  meshes.  The  first 
mesh  is  aligned  with  the  discontinuities  of  the  gradient  of  the  solution  and  the 
second  mesh  is  unstructured.  The  results  are  reported  in  Figure  2.  For  both  mesh 
types,  we  observe  that  the  approximate  L1-minimizer  is  similar  to  the  Lagrange 
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Figure  2:  Pentagon:  Aligned  unstructured  mesh  (left);  Non-aligned  unstructured 
mesh  (right). 


interpolant  of  the  exact  solution  on  the  same  mesh.  This  is  what  we  should  ex¬ 
pect  intuitively:  i.e.,  the  L1  -minimization  process  solves  the  equation  in  the  region 
where  the  solution  is  smooth  and  simply  ignores  the  PDE  in  the  regions  where  the 
gradient  of  the  exact  solution  is  discontinuous. 

3  Surface  and  image  reconstruction 

3.1  L 1  data  reconstruction 

In  geometric  modeling  and  image  reconstruction,  one  often  tries  to  extract  a  shape 
or  recover  a  piece-wise  smooth  surface  from  a  set  of  measurements.  That  is,  one 
wants  to  find  a  surface  that  satisfies  constraints  or  fits  given  data  and  is  visually 
pleasing.  The  objectives  could  vary  with  the  applications  but  the  intuitive  goal  is 
to  preserve  the  shape  of  the  object.  For  example,  one  may  want  to  reconstruct  a 
convex  body  if  the  underlying  data  comes  from  a  convex  object,  a  flat  surface  if  the 
data  is  locally  flat,  or  preserve  a  particular  structure  of  the  level  sets.  Sometimes, 
this  type  of  problem  is  solved  by  minimizing  an  L2-norm  of  the  Hessian. 

In  [15,  3],  we  have  taken  a  different  approach  that  we  think  is  well  suited  for 
man-made  surfaces,  Digital  Elevation  Models  (DEM),  and  enhancement  of  digital 
images.  Namely,  we  minimize  the  total  variation  of  the  gradient  of  a  function  con¬ 
structed  on  a  finite  element  space  satisfying  interpolatory  constraints.  Minimizing 
the  total  variation  of  the  gradient  of  a  smooth  function  amounts  to  minimizing 
the  Id-norm  of  its  second  derivatives.  The  key  observation  is  that  using  the  Ll- 
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norm  in  the  minimization  process  produces  oscillation  free  surfaces.  In  [15,  3] 
we  have  developed  a  surface  reconstruction  technique  based  on  the  minimization 
of  the  total  variation  of  the  gradient.  Convergence  of  the  method  is  established 
and  an  interior  point  algorithm  solving  the  associated  linear  programming  prob¬ 
lem  is  introduced.  The  reconstruction  algorithm  is  illustrated  on  various  test  cases 
including  natural  and  urban  terrain  data,  and  enhancement  of  low-resolution  or 
aliased  images. 


Figure  3:  Barton  Creek  data  set.  Qi  interpolant  (top  left);  L1  -reconstruction  (top 
right);  zoom  of  the  Qi  interpolant  (bottom  left);  zoom  of  the  L1 -reconstruction 
(bottom  right). 

In  Figure  3  we  reconstruct  the  elevation  map  of  an  area  near  Barton  Creek  in 
Austin,  Texas.  We  use  the  elevation  map  of  a  3km  x  3km  terrain.  The  data  comes 
from  the  Digital  Elevation  Models  (DEM)  data  files  produced  by  the  U.S.  Geolog¬ 
ical  Survey  (USGS)  (DEM  data  is  available  at  http://www.webgis.com/terraindata.html). 
The  data  set  in  sampled  on  100  x  100  uniform  Cartesian  mesh.  The  QI  interpolant 
of  the  data  and  the  L 1  cubic  reconstruction  are  shown  in  Figure  3  (two  top  panels). 

To  better  appreciate  the  quality  of  the  L 1  cubic  reconstruction  we  show  a  zoom  of 
a  small  region  in  the  two  bottom  panels  of  Figure  3. 

In  Figure  4  we  show  how  the  L 1  reconstruction  technique  can  help  to  recon¬ 
struct  blurred  or  aliased  text.  The  original  aliased  picture  is  in  the  top  panel  of 
Figure  4.  It  is  the  word  CESKOSLOVENSKO.  The  left  panel  in  the  second  and 
third  row  are  zooms  of  the  original  picture.  In  the  second  row  the  image  is  in  gray 
scale  and  in  the  third  row  we  use  only  black  and  white.  The  center  and  right  pan- 
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(a)  Original  picture. 


(b)  Zoom  of  original  picture. 


MUNI 


(c)  L1  -reconstruction. 


Figure  4:  CESKOSLOVENSKO  test  case.  Original  picture  (left  column).  Zoom 
of  original  picture  (middle  column,  gray  image  (top),  BW  thresholding  (bottom)). 
L1 -reconstruction  (right  column,  gray  image  (top),  BW  thresholding  (bottom)) 


els  are  two  different  cubic  L 1  reconstructions.  It  is  clear  that  our  technique  really 
improves  the  image  since  one  can  certainly  read  the  aliased  word  in  the  picture  in 
the  bottom  right  panel. 


(a)  Original  picture  zoomed  (b)  L1  -reconstruction 


Figure  5:  Pepper  test  case. 

We  now  show  how  our  L1  -reconstruction  method  works  on  a  pepper  image.  A 
zoom  of  original  down-sampled  image  is  shown  in  the  left  panel  in  Figure  5.  The 
result  obtained  using  the  L1  -reconstruction  is  shown  in  the  right  panel  in  Figure  5. 
The  improvement  is  clear  again. 
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3.2  Local  L1  -minimization  and  parallelization 

The  key  question  one  faces  when  using  L1  -minimization  techniques  is  to  make 
the  minimization  algorithm  fast  and  capable  of  solving  very  large  problems.  We 
have  considered  two  possible  approaches:  (i)  local  L1 -minimization,  which  allows 
parallel  implementation;  (ii)  different  implementation  method  (i.e.,  augmented 
Lagrangian  algorithm).  Both  of  these  were  implemented  for  image/surface  en¬ 
hancement  and  linear  transport  problems. 


Figure  6:  Lena  Test.  Original  image  (top  left);  down-sampled  image  (top  center); 
Standard  bi-cubic  reconstruction  (top  right);  First  Jacobi  iteration  of  local  Ll- 
reconstruction  (bottom  left);  Third  Jacobi  iteration  (bottom  center);  Global  Ll- 
reconstruction  (bottom  right);. 

In  our  previous  works  the  L1  -minimization  problems  were  formulated  glob¬ 
ally  over  the  whole  domain  of  interest,  see  [2,  3]  and  even  though  the  results 
obtained  are  excellent,  the  computational  time  are  too  large  for  large  problems. 
To  address  this  issue  we  have  developed  a  new  method  based  on  a  divide-and- 
conquer  strategy.  Instead  of  computing  the  minimizer  over  the  whole  domain  at 
once,  we  divide  the  domain  of  interest  into  sub-domains  and  compute  local  L 1 


7 


minimizers  in  each  sub-domain.  This  local  L 1  minimizing  technique  is  used  in  a 
Jacobi  iterative  scheme  to  reach  the  global  minimizer.  The  Jacobi  iterations  are 
repeated  until  a  residual  is  within  a  fixed  tolerance.  The  construction  of  the  local 
L 1  minimization  problems  is  not  trivial  and  it  is  not  evident  that  by  performing 
local  minimizations  one  can  reach  the  global  L1  -minimizer.  Our  numerical  results 
confirms  that  this  is  true  in  all  the  test  cases  that  we  have  considered  so  far. 

Figure  6  shows  results  for  the  standard  Lena  test  problem.  We  have  down- 
sampled  the  512x512  gray-scale  original  image  to  a  128  x  128  image  by  averaging 
4x4  pixel  blocks,  and  then  the  down-sampled  image  is  reconstructed  with  both 
the  local  and  the  global  L1  minimization  algorithm.  The  improvement  over  the 
standard  bi-cubic  algorithm  used  in  practice  is  clear.  These  tests  reveal  that  1) 
local  solutions  on  aggregates  are  faster  than  solving  the  global  problem;  2)  only  a 
few  (in  many  cases  just  one)  Jacobi  iterations  are  enough  to  be  close  to  the  global 
L1-minimizer  for  practical  purpose.  This  technique  can  used  for  designing  fast 
algorithms  for  surface  reconstruction  and  transport  equations  with  sparse  local 
structures/residuals.  Our  numerical  tests  show  that  the  local  algorithm  performs 
well  on  the  Digital  Elevation  Models  that  we  considered  in  [2,  3]. 


Figure  7:  Pepper  test.  Down-sampled  image  (top  left);  standard  Bi-cubic  recon¬ 
struction  (top  right);  First  Jacobi  iteration  of  local  L1 -reconstruction  (bottom  left); 
global  L1 -reconstruction  (bottom  right);. 
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The  local  L 1  minimization  algorithm  can  be  implemented  in  parallel.  In  a 
parallel  setting  the  local  sub-problems  are  solved  concurrently,  and  the  final  result 
is  communicated  between  CPUs  afterwards.  Table  1  shows  the  speed-up  obtained 
on  a  machine  with  64  processors  for  the  Lena  and  the  pepper  tests.  Note  that  these 
speed  up  figures  are  conservative,  as  we  obtain  them  by  performing  a  few  more 
Jacobi  iterations  than  one  typically  needs  in  practice. 


Problem 

Global  algorithm 

Local  parallel  algorithm 

Speed  up 

Lenna  test 

2300  s 

56  s 

41 

The  pepper  test 

1741  s 

21  s 

83 

Table  1:  CPU  time  for  the  global  and  local  parallel  algorithms  (64  CPUs). 


Another  direction  of  improvement  is  to  use  a  different  algorithm  than  the  inte¬ 
rior  point  method  for  computing  the  L 1  minimizers.  We  have  explored  an  operator 
splitting  technique  known  as  the  Augmented  Lagrangian  algorithm.  We  have  been 
able  to  show  that,  depending  on  the  accuracy  required,  the  Augmented  Lagrangian 
can  be  faster.  Table  2  compares  performance  of  the  two  algorithms  on  the  pepper 
image  problem  for  different  accuracy  requirements,  see  Figure  7. 


tolerance 

0.5 

0.2 

0.1 

0.01 

0.001 

Augmented  Lagrangian  Algorithm 

51 

185 

612 

4250 

10000+ 

Interior  Point  Method 

173 

284 

519 

1741 

8872 

Table  2:  Time  required  to  reach  the  same  error  tolerance  for  the  Interior  point 
point  and  the  Augmented  Lagrangian  algorithm.  Times  are  in  seconds. 

Our  numerical  tests  show  that  the  augmented  Lagrangian  method  converges 
relatively  fast  initially  but  then  slows  down.  The  interior  point  is  slow  but  superior 
when  it  comes  to  achieve  high  accuracy.  We  are  currently  working  on  combining 
the  fast,  but  inaccurate,  Augmented  Lagrangian  method  and  the  slow,  but  accu¬ 
rate,  Interior  Point  method  into  an  adaptive  algorithm.  The  above  results  will  be 
reported  in  a  forthcoming  paper  [14]. 


4  Entropy  viscosity 

We  report  in  this  section  on  our  findings  regarding  approximation  of  nonlinear 
conservation  equations  using  a  new  method  that  we  call  entropy  viscosity. 
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4.1  Entropy  viscosity  for  nonlinear  scalar  conservation  laws 

Nonlinear  conservation  equations  can  all  be  put  into  the  following  general  form: 

dtu(x ,  t )  +  V-/(m(x,  t))  —  0,  x  G  U,  t  >  0,  (2) 

with  u\t=o  =  u0  and  appropriate  boundary  conditions.  The  scalar  initial  boundary 
value  problem  has  a  unique  entropy  solution  which  satisfies  an  additional  set  of 
differential  inequalities 

dtE(u)  +  V-F(u)  <  0,  (3) 

for  any  pairs  E(u)  and  F(u)  such  that  E  is  convex  and  F'(u)  =  E'(u)f{u).  The 
function  E  is  called  entropy  and  F  is  the  associated  entropy  flux.  For  convex 
fluxes  (i.e.,  if  /  is  convex)  in  one  space  dimension  it  is  known  that  one  entropy 
pair,  for  example  the  one  generated  by  E{u)  =  \u2,  is  enough  to  select  the  unique 
entropy  solution.  Physical  systems  have  at  least  one  entropy  pair  and  the  entropy 
inequality  (3)  is  the  mathematical  form  of  the  second  principle  of  thermodynam¬ 
ics.  It  is  expected  that  the  auxiliary  inequality  (3)  serves  as  a  selection  criteria 
and  guarantees  convergence  of  the  numerical  approximation  to  the  correct  physi¬ 
cal  solution  of  the  nonlinear  system.  Therefore,  it  is  desirable  (and  necessary)  to 
somehow  incorporate  the  entropy  dissipation  (3)  in  a  numerical  scheme. 

A  traditional  way  of  selecting  the  entropy  solution  of  (2)  consists  of  adding 
viscous  dissipation 


dtue  +  V-f(ue )  -  V-(eVtfe)  =  0,  (4) 

where  e  >  0  and  it  can  be  shown  in  general  that  ue  — >■  u  when  f  — >  0.  The 
parameter  e  is  usually  taken  to  be  proportional  to  the  local  mesh  size  when  con¬ 
structing  numerical  approximation  of  (4)  and  this  limits  the  convergence  rate  to 
first-order  at  most.  The  use  of  artificial  viscosity  to  solve  nonlinear  conservation 
equations  has  been  pioneered  by  Neumann  and  Richtmyer  and  popularized  later 
by  Smagorinsky  for  LES  purposes  and  by  Ladyzenskaja  for  theoretical  purposes 
in  the  analysis  of  the  Navier-Stokes  equations.  The  early  versions  of  artificial  vis¬ 
cosities  being  overly  dissipative,  the  interest  for  these  technique  has  faded  over 
the  years,  especially  in  the  Discontinuous  Galerkin  Finite  and  Finite  Volume  liter¬ 
ature,  where  up-winding  and  limiters  have  been  shown  to  be  efficient  and  to  yield 
high-order  accuracy.  Up  to  a  few  exceptions  slope  limiting  is  a  one-dimensional 
concept  that  does  not  easily  generalize  to  unstructured  meshes  in  higher  dimen¬ 
sions.  Moreover,  the  theoretical  understanding  of  the  stability  and  convergence  of 
limiters  is  currently  restricted  to  uniform  grids  and  scalar  equations  in  one  space 
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dimension.  For  the  above  reasons  and,  among  other  things,  the  fact  that  artifi¬ 
cial  viscosities  are  easy  to  implement,  the  interest  for  artificial  viscosity  has  lately 
been  revived  in  the  DG  literature  and  in  the  Continuous  Galerkin  (CG)  literature 
as  well. 

We  have  introduced  a  new  technique  for  generating  high-order  numerical  ap¬ 
proximations  for  nonlinear  conservation  equations  in  [9,  5,  10]  using  continuous 
finite  elements  and  spectral  elements.  The  main  stabilization  mechanism  in  this 
method  is  a  nonlinear  dissipation  proportional  to  the  local  size  of  an  entropy  pro¬ 
duction.  For  this  reason  the  method  is  called  entropy  viscosity.  It  is  usually  argued 
in  the  literature  that  good  artificial  viscosities  can  be  computed  from  measures  of 
the  local  regularity  of  the  solution  or  from  the  local  residual  of  the  PDE.  We 
propose  to  take  a  slightly  different  route  by  considering  the  local  residual  of  an 
entropy  equation  to  construct  the  artificial  viscosity. 

e~  ch2\dtE(u)  +  V-F(u)|.  (5) 

One  immediate  consequence  of  this  choice  is  that  the  viscosity  is  proportional  to 
the  entropy  production,  which  is  known  to  be  large  in  shocks  and  to  be  zero  in 
contact  discontinuities.  As  a  result,  this  strategy  makes  an  automatic  distinction 
between  shocks  and  contact  discontinuities,  and  this  subtle  distinction  cannot  be 
made  by  any  of  the  two  classes  of  methods  mentioned  above.  We  also  think 
that  using  the  residual  of  the  conservation  equation  may  be  less  robust  than  using 
the  entropy  residual.  This  argument  is  based  on  the  observation  that  consistency 
requires  the  residual  of  the  PDE  to  converge  to  zero  in  the  distribution  sense  as  the 
mesh-size  goes  to  zero,  whereas  the  very  nature  of  entropy  implies  that  the  entropy 
residual  converges  to  a  Dirac  measure  supported  in  the  shocks.  This  implies  that 
the  entropy  residual  focuses  far  better  in  shocks  than  the  PDE  residual,  and  it  is 
in  this  sense  that  we  claim  that  the  PDE  residual  is  less  reliable  than  the  entropy 
residual.  This  can  be  better  understood  by  considering  the  simple  case  of  the  one¬ 
dimensional  Burgers  equation  dtu  +  \dxu2  =  0  with  initial  data  u(x,  0)  =  1  if 
x  <  0  and  u(x,  0)  =  0  otherwise.  The  entropy  solution  is  u(x,  t )  =  T—  H (x  —  \t) 
where  H  the  Heaviside  function.  One  observes  that  the  residual  of  the  equation 
is  zero,  whereas  the  entropy  residual  is  an  negative  measure:  \dtu 2  +  |<9«3  = 
—  j^S(x  —  \ t),  where  5  is  the  Dirac  measure.  This  effect  is  very  well  illustrated 
on  the  one-dimensional  Burgers  equation  over  the  interval  (0, 1)  with  initial  data 
sin(7r.x),  see  Figure  8.  We  show  in  this  figure  the  entropy  viscosity  computed  by 
using  either  the  quadratic  entropy  E{u)  =  \u2  (first  and  second  panel  from  the 
left)  or  E{u)  =  u  (third  and  fourth  panels).  Choosing  E{u)  =  u  corresponds  to 
using  the  residual  of  the  PDE  for  the  viscosity.  These  tests  show  that  choosing 
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Figure  8:  Entropy  viscosity  in  log  scale  for  the  Burgers  equation  at  t  =  0.25 

E{u)  =  |  u2  leads  to  better  focusing  of  the  viscosity  and  is  more  robust  with 
respect  to  the  multiplier  c  (see  definition  (5)). 

Although  no  convergence  proof  of  the  entropy  viscosity  method  has  been  pro¬ 
duced  yet,  the  method  has  been  shown  in  [9,  5,  10]  to  deliver  high-order  accuracy 
by  testing  it  on  a  large  variety  of  benchmark  problems.  We  have  recently  proved 
that  the  algorithm  with  explicit  time  stepping  is  indeed  L2-stable  for  nonlinear 
scalar  conservation  equations,  [1]. 

4.2  Compressible  fluid  dynamics 

We  have  generalized  the  entropy  viscosity  concept  to  the  compressible  Euler  equa¬ 
tions  by  using  the  physical  entropy  to  construct  the  entropy  viscosity. 

4.2.1  Continuous  finite  elements 

We  have  tested  the  method  with  continuous  finite  elements.  We  show  in  Figure  9 
the  density  field  at  t  =  2.86  and  the  viscosity  field  at  t  =  4  for  the  classical  wind 
tunnel  problem  with  a  forward  facing  step  at  Mach  3.  We  observe  that  the  viscosity 
focuses  very  well  and  there  is  almost  no  viscosity  in  the  contact  discontinuity  that 
develops  at  the  top  of  the  flow. 

We  show  in  Figure  10  the  density  field  at  t  =  0.2  for  the  double  Mach  reflec¬ 
tion  problem  at  Mach  10.  We  observe  again  that  there  is  almost  no  viscosity  in 
the  contact  discontinuity  that  develops  from  the  first  triple  point;  moreover  the  jet 
develops  the  standard  instability  at  the  bottom  of  the  domain. 
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(a)  Density  at  t  =  2.86  (b)  Entropy  viscosity  flux  at  t  =  4 

Figure  9:  Wind  tunnel  with  a  step  at  Mach  3,  Pi  approximation. 


(a)  Density  (b)  Entropy  viscosity 


Figure  10:  Double  Mach  reflection,  t  =  0.2,  Px  approximation. 

We  have  also  developed  two  different  techniques  to  enforces  boundary  condi¬ 
tions  on  curved  boundaries,  which  is  not  an  easy  task  for  nonlinear  conservation 
equations.  We  are  also  using  the  entropy  viscosity  framework  to  construct  a  goal 
oriented  refinement  strategy.  The  method  is  well  developed  by  now,  see  [10,  7,  8]. 
We  have  solved  various  benchmark  problems  with  finite  elements,  spectral  ele¬ 
ments  and  Fourier.  We  show  typical  results  in  Figure  1 1 . 


Figure  1 1 :  Supersonic  around  a  cylinder  in  a  tunnel.  Pi  approximation.  Density 
field  (left)  and  entropy  viscosity  (right). 
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4.2.2  Discontinuous  Galerkin  finite  elements 

The  entropy  viscosity  method  has  been  extended  to  the  Discontinuous  Galerkin 
setting  in  [16].  In  introducing  appropriate  numerical  fluxes  the  method  is  proved 
to  be  consistent  with  an  entropy  inequality.  The  numerical  performance  of  the 
method  is  observed  to  be  independent  of  the  chosen  polynomial  degree  of  approx¬ 
imation.  To  the  best  of  knowledge  the  extension  of  the  entropy  viscosity  method 
to  discontinuous  finite  elements  had  never  been  done  before. 


(a)  Qi  (b)  Q2  (c)  Q3 


Figure  12:  Riemann  problem  12  at  T  =  0.2  using  discontinuous  Galerkin  plus 
entropy  viscosity  stabilization.  Qi  (left),  Q2  (center),  Q3  (right),  128x128  cells. 

As  an  illustration  of  the  method  we  demonstrate  its  performance  in  Figure  12 
on  the  classical  Riemann  Problem  12  from  a  paper  by  Liska  and  Wendroff.  It 
is  a  two-dimensional  Riemann  problem  developing  complex  structures  involving 
shocks  and  contacts.  The  computational  domain  is  Cl  =  (0,  l)2.  The  heat  capacity 
ratio  is  7  =  1.4  and  the  initial  data  is 


p  =  1, 

II 

u  =  (0,0) 

0  <  x  <  1/2, 

0  <y<  1/2, 

p  =  1, 

P= 

u  =  (3/VT7, 0) 

0  <  x  <  1/2, 

1/2  <  y  <  1, 

p=  1, 

P=  !. 

u=  (0,3 /VT7), 

1/2  <  x  <  1, 

0  <y<  1/2, 

P  =  2/5, 

p  =  17/32, 

u  =  (0,0), 

1/2  <  x  <  1, 

0.5  <  y  <  1. 

Due  to  the  finite  speed  of  propagation  of  perturbations,  the  solution  of  the  problem 
in  (0,  l)2  is  identical  to  the  restriction  to  (0,  l)2  of  the  solution  to  the  Riemann 
problem  in  ]R2  up  to  time  t*  :=  2(g2+Q  6)  >  0.32,  where  s  =  ^=.  The  time 
stepping  is  done  with  RK4  with  CFL  =  0.25.  The  computations  are  done  with 
Qi,  Q2  and  Q3  discontinuous  finite  elements  on  a  grid  composed  of  16384  =  1282 


14 


quadrangular  cells.  The  total  number  of  scalar  degrees  of  freedom  for  the  Qi,  Q2 
and  Q3  approximations  are  65536,  146456,  262144,  respectively,  i.e.,  47x  (k  +  1). 

We  show  in  Figure  12  the  density  field  at  T  =  0.2  <  T*  for  the  Qi,  Q2  and  Q3 
approximations.  The  results  compare  well  with  reference  solutions.  The  shocks 
and  the  fine  structures  that  develop  behind  them  are  very  well  described.  The 
method  behaves  well  as  the  polynomial  degree  of  the  approximation  increases. 

4.3  Lagrangian  hydrodynamics 

Through  collaborations  with  colleagues  at  Lawrence  Livermore  National  Labo¬ 
ratory,  we  have  recently  started  to  extend  the  entropy  viscosity  methodology  to 
Lagrangian  hydrodynamics.  Our  first  experiments  have  shown  that  the  method 
extends  naturally  to  this  setting  without  any  particular  difficulty,  thereby  proving 
again  that  the  methodology  that  we  propose  is  very  flexible.  This  is  an  ongoing 
work,  no  formal  report  has  yet  been  written,  but  we  can  show  some  preliminary 
computations  done  by  a  student  that  has  been  supported  by  the  grant  (V.  Tomov). 


(a)  Riemann  12,  Q2,  32x32 


(b)  Sedov,  Q3,  32x32 


Figure  13: 

We  show  in  the  left  panel  of  Figure  13  the  solution  of  the  Riemann  12  problem 
at  t  =  0.2.  The  deformation  of  the  Lagrangian  mesh  is  clearly  visible.  This 
computation  has  been  done  with  Q2  continuous  finite  elements.  We  show  in  the 
right  panel  of  Figure  13  the  solution  of  the  so-called  Sedov  blast  problem.  Here 
again  the  deformation  of  the  Lagrangian  mesh  is  clearly  visible.  This  computation 
has  been  done  with  Q3  continuous  finite  elements. 
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